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SUMMARY 

The  classical  EHL  point  contact  problem  is  solved  using  a  new 
"system-approach, "  similar  to  that  introduced  by  Houpert  and  Hamrock  for  the 
line-contact  problem.  Introducing  a  body-fitted  coordinate  system,  the  troublesome 
free-boundary  is  transformed  to  a  fixed  domain.  The  Newton-Raphson  method  can  then 
be  used  to  determine  the  pressure  distribution  and  the  cavitation  boundary  subject  to 
the  Reynolds  boundary  condition.  This  method  provides  an  efficient  and  rigorous  way 
of  solving  the  EHL  point  contact  problem  with  the  aid  of  a  supercomputer  and  a 
promising  method  to  deal  with  the  transient  EHL  point  contact  problem.  A  typical 
pressure  distribution  and  film  thickness  profile  are  presented  and  the  minimum  film 
thicknesses  are  compared  with  the  solution  of  Hamrock  and  Dowscn.  The  details  of  the 
cavitation  boundaries  for  various  operating  parameters  are  discussed. 


NOMENCLATURE 


a,b  semi-major  and  semi-minor  axes  of  contact  ellipse,  m 
E  The  Young  modulus  of  elasticity,  Pa 
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equivalent  elastic  constant,  Pa  _ 
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normal  force,  N 

dimensionless  materia?,  parameter,  aE' 
dimensionless  cavitation  boundary 
cavitation  boundary  function 
dimensionless  film  thickness 
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H  .  dimensionless  minimum  film  thickness 

min 

HQ  dimensionless  reference  film  thickness 

h  film  thickness,  m 

h  .  minimum  film  thickness,  m 

nun 

hQ  reference  film  thickness,  m 

k  ellipticity  parameter 

n  normal  direction 

P  dimensionless  pressure 

PQ  Roelands  pressure-viscosity  constant,  l,96e8  Pa 

p  pressure,  Pa 

ph  Hertzian  pressure,  Pa 

u^  mean  entraining  velocity,  m/s 

x,y  cartesxan  coordinates 


INTRODUCTION 


In  the  design  of  nonconformal  contact  machine  elements,  knowledge  of  elasto- 
hydrodynamic  lubrication  (EHL)  is  needed.  Since  the  1970 's,  several  authors  have 
presented  their  results  of  the  point-contact  EHL  problem.  Among  them,  the  following 
Hamrock  and  Dowson  (H.D.)  formula  (ref.  1)  is  widely  used  in  the  design  of  many 
machine  elements: 


(H  in) 


min'H.D. 


3 . 63U°'68G°'49W  ”°'073 ( 1  -e°-68k) 
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For  an  EHL  solution,  a  nonlinear  integro-dif ferential  equation  must  be  solved, 
the  Reynolds  equation  and  the  elasticity  equation.  The  nonlinear-.  ,;ies  are  due  to: 

(1)  the  dependence  of  the  lubricant  properties,  (viscosity  and  density),  on  the 
pressure;  (2)  the  dependence  of  the  film  thickness  on  the  pressure;  and  (3)  the  free 
boundary  at  the  exit  region.  Even  for  the  hydrodynamic  lubrication,  since  the  free 
boundary  is  dependent  on  the  pressure  distribution,  the  Reynolds  equation  has  non¬ 
linear  characteristics.  It  is  well  known  to  computational  lubrication  engineers  that 
the  numerical  treatment  of  the  point-contact  EHL  problem  has  inherent  difficulties. 
One  can  see  how  difficulties  arise  upon  careful  consideration  of  the  above  mentioned 
nonlinearities.  For  example,  one  of  the  major  difficulties  is  the  piezoviscous 
effect.  At  high  loads  the  viscosity  of  the  fluid  can  vary  by  10  orders  of  magnitude 
within  the  conjunction,  which  caused  the  pressure  spikes  and  numerical  difficulties. 


Another  difficulty  associated  with  solving  the  EHL  point-contact  case  is  to 
locate  the  free  boundary  where  cavitation  occurs.  In  the  solution  presented  by  H.D. 
Christophergon' s  method  (ref.  2)  was  used  together  with  a  Gauss-Seidel  iterative 
scheme.  The  essence  of  this  method  is  to  truncate  the  negative  pressures  as  they 
occur  during  iteration  and  the  outlet  boundary  is  located  automatically.  Oh  and 
Rhode  (ref.  7)  solved  the  point  contact  EHL  problem  using  a  finite  element  method  and 
Newton's  method.  But,  it  has  been  found  that  the  nonnegativity  condition  was  needed 
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Newton's  method.  But,  it  has  been  found  that  the  nonnegativity  condition  was  needed 
to  be  checked  in  each  iteration  and  the  discrimination  between  the  continuous  film 
region  and  the  cavitated  region  was  troublesome.  Though  the  solution  can  be  obtained 
it  is  unavoidable  that  the  solution  is  dependent  upon  the  mesh  size  distribution  near 
the  boundary. 

Finally,  the  large  amount  of  computation  time  and  computer  memory  space  are 
concerns  in  this  calculation.  Toe  majority  of  CPU  time  is  devoted  to  the  calculation 
of  the  elastic  deformation.  In  general,  the  Gauss-Seidel  iterative  method  requires 
more  than  one  hundred  times  of  iterations  to  obtain  the  converged  solution.  Further¬ 
more,  to  obtain  a  solution  for  a  given  load,  one  additional  loop  is  required  to  find 
the  reference  film  thickness. 

This  all  adds  up  to  the  fact  that  it  is  very  difficult  to  achieve  a  stable 
solution  at  relatively  high  loads  and  short  CPU  times.  Recognizing  this,  Houpert  and 
Hamrock  (ref.  8)  devised  an  elegant  scheme  for  the  line  contact  case  that  enabled 
higher  load  calculations  and  saved  on  computational  time  as  well.  This  scheme  was  an 
adaptation  of  Okamura  (ref.  9)  and  became  known  as  the  “system  approach.”  Using  a 
Newton-Raphson  algorithm,  the  pressures,  the  integration  constant,  ard  the  reference 
film  thickness  are  found  simultaneously.  Here  advantage  has  been  taken  of  the  fact 
that  the  one-dimensional  Reynolds  equation  can  be  integrated  analytically  to  obtain 
dp/d x  and  in  turn  used  with  the  Reynolds'  boundary  conditions  to  locate  the  cavi¬ 
tation  boundary. 

To  the  author's  knowledge,  the  system-approach  has  not  been  successfully  applied 
to  the  point-contact  problem.  Unlike  the  line-contact  case,  the  two  dimensional 
Reynolds  equation  can  not  be  integrated  analytically.  However,  a  successful 
formulation  of  the  system-approach  can  nevertheless  be  accomplished  by  introducing  a 
body-fitted  coordinate  systen.  and  transforming  the  unknown  physical  boundary  into  a 
fixed  computations.*  boundary.  The  unknown  boundary  function  becomes  a  part  of  the 
system  matrix.  In  addition,  the  reference  film  tnickness  can  be  calculated 
simultaneously  as  was  done  in  the  line-contact  case.  This  reduces  the  number  of 
visits  to  the  elastic  deformation  subroutine  substantially.  However,  as  was  pointed 
out  by  Lubrecht  et  al.,  (re!.  3),  computer  memory  may  be  a  problem  since  the  Jacobian 
matrix  is  a  full  matrix  o\  *  „  » he  elasticity  equation.  This  problem  can  be  overcome 

by  using  the  block  tridiagonal  approximation  of  the  system  matrix.  The  matrix 
inversion  is  accomplished  by  the  Thomas  algorithm,  and  there  is  no  need  to  store  the 
whole  Jacobian  matrix.  Furthermore,  the  force-balance  loop  can  be  obviated  by 
including  it  in  the  system  equations  and  solving  simultaneously. 

In  this  paper,  the  classical  EHL  point-contact  problem  is  revisited  with  a  new 
formulation:  a  free  boundary  value  problem  using  the  system-approach  described 

above.  The  minimum  film  thicknesses  are  compared  with  equation  (1>  and  the  details 
of  the  free  boundary  are  discussed. 


2.  ANALYTIC  FORMULATION 
2.1  Contact  Geometry 

Figure  1  shows  the  physical  model  of  the  elliptical  contact,  where  the  x-axis 
represents  the  rolling  or  sliding  direction,  and  yg  are  the  inlet  boundaries 
and  x  =  g(y)  is  the  outlet  or  cavitation  boundary.  The  ellipticity  parameter  is 
expressed  in  terms  of  the  curvature  difference  (T),  the  elliptic  integral  of  the 
first  (F),  and  the  second  kind  (S)  as 
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.  flE? 


where 


S(1  +  T) 


S(1  -  T) 


F  -  J0  1  -  1  -  _  sin'<p\  dp 


fo  *|l  '  1  ~  ~  sin^[  d<f> 


T  Jl-1 

Rx  Ry 


Defining  r  as  R/R  ,  equation  (2)  can  be  rewritten  as 


k  *  |(r  +  1)  -  -  r 

l  s 

Therefore,  given  r,  the  ellipticity  parameter  can  be  calculated  iteratively. 


2.2  Governing  Equations 

Assuming  isothermal  conditions  and  that  the  lubricant  xs  Newtonian,  the  steady- 
state  Reynolds  equation  for  the  point  contact  problem  is 


9  ph3  9p  ^  9  /Jh3  9p  _  5(/>h) 

9x[  p  9x  9y  k  p  9y;  Um  9x 


and,  using  the  parabolic  approximation  for  the  geometry,  the  film  thickness  is 
expressed  by 


2  2 


»<x,y>  -  h„  ♦  JL  .  JL  *  _L  f  f  .  PU'.Y'Idx'dy- 

2R  2R  ms.'  J  Jn  I - 


2RX  2Ry  7T  E' 


/(x  -  x' )2  + 


(y  -  y'  )* 


The  applied  normal  force  may  be  balanced  by  the  generated  hydrodynamic  pressure 
distribution, 


f  “  J  Jn  p<x,y>dx  dy 
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Applying  the  Reynolds  boundary  condition  and  symmetry  condition  at  the  x-axis, 
the  boundary  conditions  are: 


H 

O 

at 

x  a  XA;  0  <,  y  i  y0. 

p  *  0 

at 

xA  a  $  x  <  g(yB) ;  y  = 

0p 

5S'° 

at 

x  «  g(y) ;  0  i  y  £  yB, 

!?-o 

at 

xA$  x  i  g  ( 0 ) ;  y  «  0  - 

9y 

The  viscosity-pressure  relation  is  modeled  by  the  Roelands  (11)  equation,  i.e., 

Poaf  pV 

/l  *  H0  ex p  —  1  +  —  -1  k 

2  P0 


a  *  _(ln  H0  +  9.67)  ( 

Po 

and,  the  Dowson-Higginson  relation  (12)  is  used  for  the  density-pressure  relation, 

0.59X109  +  1.34p 

P  =  Pq  - -Pin  Pa  (1 

0.59xl09  +  p 


Letting 


2.3  The  Dimensionless  Equations 


p  x  y  -  g 

P  *  _ ,  H  *  _ ,  X*_,Y=*_.,G*_ 

P*.  v.2  b  a  b 


3f  —  U  -  p 

a  -aph'  ph  -  r— /»  *  P  =  T- 


2^ab  /»0  />„ 


the  equations  (4)  to  (6)  become 


9  pH3  9p  1  9  pH2  9p  ,  9 

'^(T^J  ax  ^ 


H  ( X ,  Y )  «  H0  +  1  (x2  +  c^2) 


P(X' , Y '  )dX'  dY' 

(X  -  X'  )2  +  k2(Y  -  Y'  )2 


ft'  ,  27T 

J  JnP(X'Y)dX  dY  -  _ 


12U.u_R 

“  0  m  x 


Ct  -  C2 


k2  MPhRx 

- ,  C-  a  - 


where,  X 


E'  b 


The  dimensionless  parameters  used  here  can  be  related  to  those  used  by  Hamrock  and 
Oowson  as  follows, 


87TUk 


W 


3GWk 

a 


2ir 

kTT  ( 1  +  r ) 
4Sr 


(14) 


a 
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6Wk  S  r 


when  K 


7T  1  +  r 

1  (circular  contact),  c x  and  c2  are  1 


1/3 


2 ♦ 4  Coordinate  Transformation 

Introducing  the  body-fitted  coordinate  system  described  in  reference  13, 

YB<X  -  XA> 


t 


G  -  X. 


V  «  Y 


the  following  equations  are  obtained: 
The  Reynolds  equation 

1^(1-, G,H0) 


YB  9 

9p 

a  _ 

l  9 
+  _  _ 

/  \ 

9p 

.  e 9 

/ 

G'  9p 

rr 

(6  -  x.,!  9e 

l  Hi 

k2  ^ 

{  9j?J 

k’  dri 

G  -Xa9^ 

V 

9 


k2(G  -  X.) 


9p 
q  *- 


driJ 


e(G')z  _9 

2  9f 


k  (G  -  X.) 


,  dp 

lq<  5|J 


(15) 


(16) 


.  Yg  9  _ 

-  X  -  —  (pH)  -  0 


G  -  Xa  9£ 


where,  represents  dG/dY. 

The  film  thickness  equation 


H It'D  -  H0  ♦  _ 


C(G  ~XA) 


+  x. 


CIV4 


2c2 

+  —  D(P,G) 

n2 


(17) 


The  D  represents  an  integral  operator  which  calculates  the  elastic  deformation  of 
two  solids  in  contact  resulting  from  the  pressure  distribution  in  the  fluid  film 
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region  <Q).  In  this  paper,  the  technique  presented  by  Chang  (ref.  14)  is  used.  This 
method  provides  and  efficient  way  of  evaluating  D  without  lengthy  and  complex 
mathematical  expressions.  Since  the  coordinate  transformation  can  easily  be 
implemented  to  this  method,  the  details  of  the  algorithm  are  not  presented  here. 


The  force  balance  equation  is 

l2(p,g,  .2J:°  jv0b  p(M> 


°  '  XA  ,  27T 

-  d£  d77  -  _  «  0 

3 


(18) 


In  the  above  formulation,  hl  is  the  nonlinear  partial  differential  operator  and  L2 
is  integral  operator. 


3.  NUMERICAL  METHOD 
3.1  Spatial  Discretization 

To  provide  a  small  mesh  size  near  the  pressure  spike  region,  an  interior 
stretching  function  (ref.  15)  is  adopted  along  the  £-axis.  The  finite  difference 
representation  of  the  transformed  Reynolds  equation  is  provided  in  the  appendix. 


3.2  Newton's  Method 


The  system  equations  are 

L^( P, G, HQ)  -  0 

L2(P,G)  -  0  <19 

In  reference  13,  for  hydrodynamic  case,  has  been  solved  using  the  Thomas 
algorithm  to  find  P  and  G,  and  L2  was  used  to  determine  Ho  by  the  force 
balance  loop.  For  the  EHL  case,  the  same  method  can  be  used.  However,  if  L2  can 
be  put  in  the  system  equation  and  be  solved  simultaneously  without  creating  a 
computer  memory  problem,  the  computation  will  be  greatly  reduced. 


The  system  equation  for  Newton's  method  can  be  written  as 


[A)  {B} 

'W 

at  — 

>4 

,  (<f  o, 

W 

L2 

/ 

(u)ntl  -  (u}n  ♦  {$u} 

n+1  n  c 

H0  *  H0  +  6H0 


(20) 

(21) 

(22) 


where 

{u}  -  |pr  J  Gj),  1-2,  NX  -  1,  J  -  1,  NJ-1 
and  [A],  {B},  {C>T  are  the  elements  in  the  Jacobian  matrix. 

In  EHL,  due  to  the  integral  operator  D,  (A]  is  full  matrix.  3ut,  since  large 
amounts  of  storage  and  computational  time  are  required  to  solve  it,  the  block 
tridiagonal  matrix  approximation  can  be  used.  Each  block  matrix  is  a  full  matrix 
which  is  different  from  the  one-sided  arrow  shape  matrix  that  resulted  from  the 
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hydrodynamic  case  where  the  elasticity  equation  is  not  needed.  The  unknown  G 
matches  the  residual  function  at  the  cavitation  boundary  where  the  pressures  are 
known  from  the  Reynolds  boundary  condition. 


Equation  (20)  can  be  rearranged  as 


[A)N  *  {b}$h0  -  -(lx) 

(23) 

{c)T(5u}  -  l2 

(24) 

From  equation  (23) 

{6u}  -  -[aj-,(l,}  -  [A)-1(b}<5h0 

(25) 

Then  using  equations  (24) 

and  (25) 

l.-tfw-fr,} 

(26) 

And  {5u>  can  be  calculated 

{cftAr’jB} 

using  equations  (25)  and  (26). 

In  equation  (26),  ( A  J _1 <  B  > 

and  (A]“1{Ll>  are  obtained  by  the  Thomas  algorithm  and  then  stored  temporarily  and 
used  in  equation  (25). 


The  convergence  criteria  are 


(1)  Pressure 


(2)  Cavitation  boundary 


EE 

i  j 


i  +i 

0 

'  I,  J 


E  E 

I  J 


<  5.0x10 


<  5*0x10 


(3)  Reference  film  thickness 


n 


5.0xl0*3 


4.  RESULTS  AND  DISCUSSION 

The  dimensionless  material  parameter  used  in  this  analysis  is  G  =  3488  in  which 
z  *  0.55,  H0  =  0.018  Ns/m2,  U  =  0.3  and  E'  =  2.19xl0*u  N/m2.  In  figures  2  and  3, 
the  pressure  distribution  and  film  thickness  profile  for  the  circular  contact  is 
presented.  The  inlet  boundary  used  for  this  analysis  is  defined  as  =  -4.0  and 
Yfi  =  2.0.  The  maximum  pressure  is  1.33  times  the  maximum  Hertzian  pressure  or 
515  MPa  which  occurs  on  the  x-axis.  The  dimensionless  minimum  film  thickness  is  0.27 
and  it  occurs  at  the  side-lobes,  X  =  0.49  and  Y  =  0.6. 
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The  majority  of  the  computation  time  is  used  for  the  calculation  of  the  elastic 
deformations  and  the  differentiations  of  the  residual  functions  with  respect  to  the 
cavitation  boundary  function  since  the  integral  operator  is  a  function  of  it.  It 
takes  about  20  sec  on  the  CRAY-XMP  at  NASA  Lewi9  Research  Center  for  1  Newton 
iteration  with  3060  nodal  points  of  the  whole  domain,  and,  in  general,  the  converged 
s.1  .ution  can  be  obtained  within  3  iterations  as  long  as  the  initial  guess  Is  within 
the  sphere  of  attraction.  It  was  .sported  (ref.  3)  that  using  the  multigrid 
interative  method  it  took  2  hr  of  C.rU  time  on  a  VAX11/750  with  2937  nodal  points. 

Since  a  different  computer  was  used,  c  direct  comparison  is  difficult.  The  current 
method  is  quite  fast  partly  because  the  direct  matrix  inversion  of  the  block  matrices 
is  vertorizable  which  makes  it  well  suited  to  the  supercomputing.  Also  because  the 
amount  of  visits  to  the  elasticity  s  'broutine  is  small  and  there  is  no  need  of  a 
force  balance  loop.  When  the  current  work  is  used  for  transient  calculations,  the 
previous  solution  is  used  as  a  guess  to  the  next  time  step  and  it  accelerate  the 
solution  process,  but  this  is  not  true  for  the  iterative  method.  This  fact  supports 
the  current  work  as  a  good  candidate  for  transient  EHL  point  contact  computation. 

The  calculated  minimum  film  thickness  in  this  investigation  for  various  operating 
parameters  are  provided  in  table  I  along  with  those  ootained  from  the  H.D.  Formula, 
equation  1.  In  general,  the  results  from  this  analysis  were  higher  than  those 
predicted  by  the  H.D.  for  the  circular  contact  case.  However,  for  the  elliptical 
contact  our  results  were  lower.  But  the  differences  do  not  exceed  10  percent. 

Figure  4  shows  another  pressure  distribution  for  circular  contact  where  a  very 
steep  pressure  spike  occurs.  The  operating  condition  is  W  =  9.154xl0’8, 

U  =  1.62xl0~u,  or  5  *  5.723  and  X  =  0.862.  The  maximum  pressure  is  2.89  times  the 
maximum  Hertzian  pressure  or  1.04  GPa.  To  the  authors'  experience,  the  solution  is 
so  unstable  beyond  this  operation  range  that  the  convergence  usually  fails.  When 
U  *  6.432xl0“12,  the  maximum  possible  W  is  2.367xlO~7  for  circular  contact,  or  3  ■ 
7.862  and  X  *  0.0964.  According  to  numerous  computations,  it  is  found  that  the 
value  of  3  dictates  the  numerical  stability  of  current  method.  The  numerical 
stability  may  be  enhanced  by  reducing  the  step  sizes  near  the  pressure  spike  region. 
But  it  should  be  noticed  that  the  Roelands  viscosity-pressure  relation  is  known  to  be 
valid  up  to  1  GPa  or  lower.  At  such  a^high  pressure  the  lubricant  behaves  as  a 
solid-like  material  and  becomes  non-Newtonian.  Also,  recently,  it  was  observed  that 
slippage  of  the  lubricant  occurs  at  or  very  near  the  surface  (ref.  16).  Thus  it  is 
believed  that  the  modification  of  the  classical  Reynolds  equation  including  the  non- 
Newtonian  effect  and  a  more  realistic  pressure-viscosity  relation  including  the 
thermal  effect  are  needed  to  investigate  the  lubrication  performance  for  the  high 
load  and  high  speed  cases. 

Figure  5  depicts  the  calculated  free  boundaries  of  the  circular  contacts  for 
various  operating  parameters.  The  x-axis  is  somewhat  stretched  to  exaggerate  the 
differences.  The  dotted  line  represents  the  Hertzian  circle.  The  general  trend  is 
that,  as  expected,  the  low  speed  condition  results  in  a  curve  that  conforms  more  to 
Hertzian  (dry)  contact  circle.  Comparing  curves  1  and  2,  the  boundary  on  the  X-axis 
stretches  more  outside  for  the  higher  load  but  elsewhere  it  is  closer  to  the  Hertzian 
contact  circle.  Comparing  curves  3  and  4,  increasing  the  speed  parameter  leads  to  a 
thicker  film  and  tends  to  straighten  out  the  boundary. 


5.  CONCLUSIONS 

The  classical  EHL  point  contact  problem  is  solved  using  a  new  "system-approach,” 
similar  to  that  introduced  by  Houpert  and  Hamrock  for  the  line-contact  problem.  This 
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requires  inverting  a  system-matrix  (i.e.,  the  Jacobian)  which  via  a  body-fitted 
coordinate  transformation  includes  boundary  conditions  at  the  free  boundary. 

Further,  a  force-balance  loop  is  avoided.  Using  a  Newton-Raphson  algorithm,  the 
pressures,  the  cavitation  boundary  curve,  and  the  reference  film  thickness  are  found 
simultaneously.  The  method  is  computationally  fast  and  has  no  problem  with  locating 
the  cavitation  boundary.  This  study  revealed  that 

1.  The  minimum  film  thickness  obtained  in  thit  study  were  all  within  10  percent 
of  the  predictions  using  the  H.D.  Formula. 

2.  The  algorithm  is  well  suited  to  performing  transient  EHL  calculations  using 
the  supercomputer  and  the  solution  at  each  time  step  accelerates  the  succeeding 
solution. 

3.  Numerical  instabilities  were  encountered  when  the  value  of  5,  that  is,  W  or 
C  is  high.  To  obtain  a  more  stable  solution,  it  is  believed  that  the  Reynolds 
equation  should  be  modified  to  include  the  non-Newtonian  effect  and  a  more  realistic 
pressure-viscosity  relation  for  high  pressure. 

4.  The  calculated  cavitation  boundary  is  near  the  Hertzian  contact  circle  but 
deviates  it  for  high  speed. 
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where 


APPENDIX  -  DISCRETIZATION  ON  THE  TRANSFORMED  REYNOLDS  EQUATION  FOR  k  =  1. 

*  (Ri  +  ^i*i/2R7 ) — - ~  pi,j)  ~  (R 1  +  Ci-i/:R7)(3i-i/2,j<pi,j  “  pr-i,j) 
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qI,J 
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DP. 
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